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1.  Introduction 


In  support  of  the  larger  goal  to  provide  a  plasma  engineering  capability  to  the  spacecraft 
community,  the  objective  of  the  Spacecraft  Charging  Modeling — Nascap-2k  contract  is  to 
develop,  incorporate,  test,  and  validate  new  algorithms  for  the  three-dimensional  plasma- 
environment  spacecraft  interactions  computational  tool,  Nascap-2k.  Nascap-2k  is  being  modified 
to  extend  the  range  of  plasma  physics  phenomena  that  the  code  can  simulate,  make  the  advanced 
code  capabilities  more  accessible  to  users,  and  improve  and  maintain  both  the  graphical  and  non- 
graphical  interfaces  to  the  code.  The  upgraded  code  is  being  used  to  simulate  problems  of 
interest  to  AFRL,  including  support  of  the  DSX  (Demonstration  and  Science  experiments) 
spacecraft  and  program. 

1.1.  Third  Year  Progress 

Details  on  progress  during  the  period  from  September  21,  2013  through  September  19,  2014  are 
given  in  the  following  sections.  Key  accomplishments  included  in  Version  4.2  are  the  following: 

•  Nascap-2k  has  been  made  significantly  more  user  friendly.  The  scripting  was  revised  to 
make  the  commands  more  intuitive.  The  interface  was  modified  to  be  more  automated 
and  intuitive. 

•  The  Nascap-2k  charging  capability  was  generalized  to  address  a  wider  array  of  charging 
problems,  including  charging  in  dense  plasma. 

•  Nascap-2k  4.2  was  prepared  for  release,  including  updating  the  documentation. 

1.1.1.  Maintenance  and  Support 

The  Nascap-2k  scripting  for  Version  4.2  was  revised  to  make  the  functionality  of  the  commands 
more  intuitive: 

•  A  new  primary  script  command  Read  Object  was  added  to  replace  the  “ReadObject” 
subcommand  of  the  Charge  Surfaces  script  command.  It  also  initializes  surface 
potentials. 

•  The  subcommands  of  the  Charge  Surfaces  command  that  were  used  to  perform 
initialization  and  database  access  operations  were  eliminated.  These  functions  are  now 
performed  automatically  as  needed. 

•  The  subcommands  used  to  set  calculational  parameters  were  combined  into  a  single  “Set 
Parameters”  subcommand. 

•  Old  scripts  still  run  correctly. 

•  Additional  discussion  of  these  revisions  appears  in  Section  2. 

The  Nascap-2k  charging  capability  for  Version  4.2  was  generalized: 

•  The  graphical  user  interface  and  the  default  scripts  were  revised  to  allow  for  a  wider 
variety  of  charging  problems  to  be  easily  specified.  The  underlying  coding  was 
streamlined  at  the  same  time. 
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•  The  ability  to  compute  and  use  secondary  electron  currents  in  all  charging  calculations 
with  tracked  currents  was  added.  The  yield  is  computed  using  analytic  currents  with  the 
correct  ion  masses. 

•  The  portions  of  the  code  that  compute  the  current  to  a  surface  and  its  derivative  with 
respect  to  potential  were  revised  to  be  more  transparent  and  easier  to  maintain. 

Comments  were  added,  the  code  reorganized,  and  duplicate  coding  removed. 

•  Additional  discussion  of  the  generalization  of  the  Nascap-2k  charging  capability  appears 
in  Section  3. 

The  Version  4.2  interface  was  modified  to  make  the  code  easier  to  use. 

•  The  ability  to  display  the  differential  potential  on  the  Results  and  Results  3D  tabs  was 
added. 

•  Earlier  a  temporary  button  had  been  added  to  the  Nascap-2k  Edit  Materials  dialog  box 
that  allowed  the  user  to  reset  the  material  properties  to  their  default  values.  The  label  on 
this  button  was  changed  to  “Reset  to  Default  Properties”  and  its  functionality  was 
expanded  to  also  immediately  redisplay  the  dialog  box  with  the  revised  values. 

•  The  lists  of  quantities  available  on  the  pull  down  menus  on  the  Results  and  Results  3D 
tabs  were  made  dynamic.  The  items  in  the  pull  down  menus  are  now  only  those  quantities 
available  in  the  database  for  the  presently  displayed  time  step. 

•  The  names  of  the  more  obscure  spatial  quantities  available  for  plotting  were  revisited  and 
the  names  of  several  quantities  on  the  pull  down  menu  on  the  Results  3D  tab  were 
revised. 

•  Object  Toolkit,  GridTool,  and  the  Nascap-2k  graphical  user  interface  were  modified  to  be 
more  tightly  coupled. 

o  GridTool  and  Object  Toolkit  now  automatically  load  the  object  when  launched 
from  the  Nascap-2k  graphical  user  interface. 

o  Nascap-2k  now  automatically  loads  the  revised  object  when  returning  from 
Object  Toolkit. 

o  When  called  from  Nascap-2k,  GridTool  now  creates  the  grid  file  with  the  right 
name  in  the  right  place  by  default. 

•  When  an  object  is  read  into  Nascap-2k,  all  of  the  pull  down  lists  of  materials  and 
conductors  are  immediately  updated. 

•  The  number  of  iterations  and  the  number  of  timesteps  are  now  the  same  variable. 

•  The  underlying  code  and  the  user  interface  were  modified  so  that  the  units  in  which 
charge  density  is  displayed  on  the  Nascap-2k  Results  3D  tab  is  now  number  per  cubic 
meter. 

•  The  Nascap-2k  user  interface  was  modified  so  that  potentials  specified  on  the  Applied 
Potentials  tab  propagate  to  the  first  instance  in  the  script  and  the  reverse. 

•  The  ability  to  display  the  time  history  of  current  components  was  added  to  Results  tab 
coding. 

Additional  improvements  to  the  default  environments,  the  algorithms,  and  the  output  files  made 
during  this  period  include  the  following: 
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•  A  new  predefined  geosynchronous  environment  based  on  an  environment  measured  by 
SCATHA  was  added. 

•  The  predefined  Freja  auroral  environments  were  modified  to  exclude  the  unphysical  low 
energy  power  law  electrons. 

•  In  Object  Toolkit,  nodes  of  the  same  component  that  are  within  a  specified  distance  of 
each  other  are  automatically  combined.  This  distance  was  changed  from  a  fixed  value  of 
one  millimeter  (0.001  m)  to  a  user  specifiable  parameter. 

•  The  maximum  length  of  a  database  name  was  determined  to  be  70  characters. 

•  The  charging  code  was  modified  so  that  charging  calculations  will  execute  without  an 
environment  and  even  without  an  object. 

•  The  script  interpretation  code  was  modified  so  that  multiple  “TimeParams”  folders  in  a 
“DoTimeSteps”  subcommand  can  be  used  to  specify  a  sequence  of  steps  with  different 
timestep  length  parameters  to  be  used  in  a  single  charging  calculation. 

•  The  values  printed  in  Tracker  text  output  files  in  reverse  trajectory  calculations  were 
modified  so  that  they  have  physical  meaning.  Macroparticles  carry  current  without  the 
Maxwellian  exponential  factor.  The  values  printed  as  the  current  to  object,  lost  current, 
and  collected  current  now  include  an  exponential  factor  appropriate  to  the  final  velocity 
of  the  macroparticle.  This  means  that  the  “lost  current”  is  the  current  actually  collected  by 
the  object  and  the  other  values  are  related  to  the  current. 

•  The  particle  specification  code  was  modified  so  that  the  emitter  current  is  assigned  to  the 
surface  element  which  is  emitting  the  current. 

•  The  particle  specification  code  was  modified  to  allow  for  injected  current  from  the 
computational  grid  boundary  in  particle-in-cell  (PIC)  calculations  when  the  electric  field 
is  outward. 

•  The  original  implementation  of  the  zero  total  current  algorithm  is  only  effective  if  all  the 
currents  are  analytic.  The  implementation  was  revised  so  that  the  algorithm  can  be  used 
in  charging  calculations  with  tracked  currents. 

Nascap-2k  4.2  is  being  prepared  for  release: 

•  The  new  Nascap-2k  functionality  was  fully  tested.  A  number  of  small  issues  were 
identified  and  the  code  was  modified  as  appropriate.  A  discussion  of  this  testing  is 
included  in  Section  4.  Example  problems  using  the  new  emitter  and  detector  capabilities 
are  included  in  Sections  4.8  and  4.9  below. 

•  The  Nascap-2k  documentation,  including  the  examples  in  Nascap-2k  User  s  Manual  and 
the  on-line  documentation,  was  updated.  It  is  now  undergoing  final  revision. 

•  The  JavaChart  source  code  (used  to  create  and  display  the  graph  on  the  Results  tab)  was 
added  to  the  Nascap-2k user  interface  Eclipse  project  and  archived  in  the  Source  Safe 
database.  This  eliminates  the  extra  javachart. jar  file  included  in  the  distribution  and  will 
insure  that  this  code  will  remain  compatible  with  future  revisions  of  Java. 

•  The  Nascap-2k  4.2  Beta  code  was  sent  to  AFRL  and  three  additional  beta  testers:  Dr. 
Michelle  Donegan  of  Applied  Physics  Laboratory  Johns  Hopkins  University,  Wousik 
Kim  of  the  Jet  Propulsion  Laboratory,  and  Nikki  Noushkam  of  Orbital  Sciences 
Corporation. 

Development  of  Nascap-2k  4.3  has  begun  with  the  following  two  activities: 
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•  Under  separate  funding,  the  potential  solution  module  was  parallelized.  As  part  of  that 
effort,  the  code  was  modified  to  dramatically  reduce  reading  and  writing  to  the  database. 
This  change  speeds  up  the  code  from  negligible  to  50%  faster  depending  on  the  number 
of  special  elements,  the  number  of  grids,  and  the  complexity  of  the  grid  structure. 

•  We  moved  the  Nascap-2k  source  code  to  the  latest  development  environment  (Microsoft 
Visual  Studio  2013)  and  compilers  (Microsoft  Visual  C++  2013  and  Intel  Visual  Fortran 
14).  We  successfully  ran  the  standard  test  cases.  With  the  new  compilers,  the  code  runs 
about  10%  faster. 

1.1.2.  Charging  Index 

Resources  regarding  internal  charging  environments  and  internal  charging  computation  were 
identified  and  summarized.  This  memo  is  included  as  Section  5  below. 

Resources  regarding  auroral  charging  environments  and  auroral  charging  computation  were 
identified  and  summarized.  This  memo  is  included  as  Section  6  below. 

We  began  an  investigation  of  DMSP  charging  events.  We  reviewed  the  literature  for  interesting 
charging  events.  We  requested  and  received  environment  and  operations  data  for  two  events.  We 
examined  the  time  history  for  these  events. 

1. 1.2.1.  SEE  Spacecraft  Charging  Handbook 

In  order  to  support  its  use  in  the  AFRL  Charging  Index  study,  we  modified  the  NASA  code,  SEE 
Spacecraft  Charging  Handbook,  so  the  code  does  not  fail  on  browser  version  verification.  Layout 
issues  on  the  3-D  charging  pages  remain.  The  browser-based  code  appears  to  work  similarly  in 
Internet  Explorer,  FireFox,  and  Opera,  but  fails  in  Safari  and  Chrome.  Note  that  the  SEE 
Spacecraft  Charging  Handbook  is  a  NASA  code  that  is  distributed  by  NASA. 

Due  to  changes  in  Java  security,  the  code  requires  the  Java  Control  Panel  security  level  be  set  to 
“Medium.”  However,  starting  with  Java  version  8u20  the  “Medium”  setting  is  no  longer 
available.  An  alternative  workaround  (for  today,  at  least)  is  to  add  the  FILE  protocol  to  the 
“Exception  Site  List,”  as  shown  in  Figure  1 . 
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Slava  Control  Panel 

leral  I  Update  j  Java  Security  J  Advanced 
J]  Enable  Java  content  in  the  browser 


Security  level  for  applications  not  on  the  Exception  Site  list 

Very  High  -  Only  Java  applications  identified  by  a  certificate  from  a  trusted  authority  are 
allowed  to  run,  and  only  if  the  certificate  can  be  verified  as  not  revoked. 


o  High  -  Java  applications  identified  by  a  certificate  from  a  trusted  authority  are  allowed  to 
run,  even  if  the  revocation  status  of  the  certificate  cannot  be  verified. 


Exception  Site  List 

Applications  launched  from  the  sites  listed  below  will  be  allowed  to  run  after  the  appropriate  security 
prompts. _ 


file:/ 


Edit  Site  List.. 


Restore  Security  Prompts 


Manage  Certificates.. 


£pplv 


|&|  Exception  Site  List 


Applications  bundled  from  the  sites  listed  below  will  be  allowed  to  run  after  the  appropriate  security  prompts. 


FILE  and  HTTP  protocols  are  considered  a  security  risk.  We  recommend  using  HTTPS  sites  where  available. 

(  <*  1  I.  Cancel  J 


Figure  1.  Adding  the  FILE  protocol  to  the  “Exception  Site  List.” 

1. 1.2.2.  Charging  in  a  LEO  Plasma 

We  explored  the  limitations  and  overlap  between  the  various  methods  of  modeling  the  charging 
of  spacecraft  with  a  low-Earth-orbit  plasma  within  Nascap-2k.  This  work  resulted  in  a  paper 
presented  at  the  Spacecraft  Charging  Technology  Conference  in  Pasadena,  California.  The  paper 
included  in  the  proceedings  is  reproduced  below  in  Section  7. 

1.2.  Reports  and  Meetings 

We  participated  in  an  in-person  contract  review  at  Kirtland  AFB  on  May  14,  2014. 

We  attended  the  Spacecraft  Charging  Technology  Conference  in  Pasadena,  California,  June  23- 
27,  2014.  We  made  an  oral  presentation  titled  Semi-analytic  and  PIC  (Particle- in-cell)  methods 
for  Quantifying  Charging  in  Dense,  Cold  Plasma.  The  paper  is  included  in  the  conference 
proceedings  and  is  reproduced  in  Section  7  below. 

These  presentations  are  included  in  the  quarterly  reports  for  the  relevant  periods. 

1.3.  Scientific  Reports  and  Journal  Articles 

An  oral  presentation  titled  “Semi-analytic  and  PIC  (Particle-in-cell)  Methods  for  Quantifying 
Charging  in  Dense,  Cold  Plasma”  was  prepared  and  presented  at  the  13th  Spacecraft  Charging 
Technology  Conference  in  Pasadena,  CA  in  June  2014.  The  paper  is  included  in  the  conference 
proceedings  and  is  reproduced  in  Section  7  below. 
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1.4.  Personnel 


Dr.  Victoria  A.  Davis  is  the  project  manager  and  Dr.  Myron  J.  Mandell  is  the  principal 
investigator.  The  following  Leidos  staff  members  have  contributed  to  the  work  reported  here. 

•  Dr.  Victoria  A.  Davis 

•  Dr.  Myron  J.  Mandell 

•  Ms.  Barbara  M.  Gardner 

•  Ms.  Lisa  Lattarulo 

2.  REVISIONS  TO  NASCAP-2K  SCRIPTING 

The  Nascap-2k  scripting  was  revised  to  eliminate  obtuse  commands.  The  changes  were  made 
backward  compatible,  so  that  old  scripts  still  run  correctly. 

The  object  is  now  read  from  the  XML  file,  placed  in  the  database,  and  the  surface  potentials  set 
to  their  initial  values  (from  the  XML  file)  by  a  new  primary  script  command  Read  Object  rather 
than  by  a  subcommand  of  the  Charge  Surfaces  script  command.  This  dramatically  reduces  the 
complexity  of  scripts  for  calculations  with  fixed  surface  potentials. 

The  obscure  subcommands  of  the  Charge  Surfaces  command,  “Calculate_Matrices”, 
“Definelnsulators”,  “InitializeCalculation”,  “OpenDatabase”,  “PrepareChargeMatrix”, 
“Read_Matrices”,  “ReadObject”,  “ReadPotentials”,  and  “WritePotentials”,  were  eliminated. 
Initialization  and  database  access  operations  are  now  performed  automatically  as  needed. 

The  subcommands  used  to  set  calculational  parameters,  “Compute  Transverse  Currents”, 
“FieldsFromFile”,  “SpaceChargeLimitedPhotoemission”,  “ZeroCurDerivAlgorithm”,  and 
“ZeroTotCurAlgorithm”  were  combined  into  a  single  “Set  Parameters”  subcommand. 

Scripts  for  common  calculations  are  now  simpler.  Figure  2  through  Figure  4  show  the  more 
compact  scripts  for  some  common  calculations. 
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GEO  Charging 


Old  Script 


New  Script 


Scnpt - 

i  r 

_ Value 

A  xrnlns 

A  Prefix  G 

o-  C  ReadObjed 
C  OpenDataoase 
C  imtiaUzeCalculatJon 
o-  C  Setlllummalion 
©-  C  SetEnvironment 
o-  C  SetlnitialCondudorPo 
C  SednitialPotentials 
C  PrepareChargeMatrix 
9  C  DoTimeSteps 

A  CloseFiles  Fi 

®-  TimeParams 
G  SetEnvironment 
9  C  DoTimeSteps 

A  CloseFiles  Fi 

~J  TimeParams 
C  SetEnvironment 
?  G  DoTimeSteps 

A  CloseFiles  Ti 

o- !<j  TimeParams 


Script - 

_ _ _ _ ( _ Value 

c  Read_Objed 

9  G  Charge.Surfaces 
o~  C  Setlllummalion 
©-  C  SetEnvironment 
C  SetlnitialCondudorPo 
9  C  DoTimeSteps 
®-  3  TimeParams 
o-  C  SetEnvironment 
9  C  DoTimeSteps 
►  TimeParams 
o-  C  SetEnvironment 
9  C  DoTimeSteps 

TimeParams 


Figure  2.  Comparison  of  Scripts  Used  to  Compute  Charging  in  a  Geosynchronous  Substorm  using  the  Old 
and  New  Schemes 


Potentials 
Old  Script 


in  Space 
New  Script 


Script - 

Value 

?  C  Charge_Surfaces 

A  xrnlns 

A  Prefix  Bipolar 

C  ReadObjed 
C  OpenDatabase 
©-  C  SetCondudorBias 
o-  C  SetCondudorBias 
C  SetCondudorBias 
C  SetCondudorBias 
o-  C  SetCondudorBias 
c  SetCondudorBias 
*■  C  SetCondudorBias 
*■  C  SetlnitialCondudorPo 
C  SetlnitialPotentials 
C  WritePotentials 
o-  C  Embed_Objed_in_Grid 
o-  C  PotentialsJn_Space 
«-  C  Create_Particles 
C  TracK_Paitides 


Scnpt - 

_ _ j  Value 

9  C  Read_Objed 

A  Filename  BipolarObjed  xml 

C  Embed_Objed_in_Gnd 
®-  C  PotentialsJn_Space 
o-  C  Create_Partides 
o-  C  Track_Partides 


Figure  3.  Comparison  of  Scripts  Used  to  Compute  Potentials  in  Space  Given  Specified  Surface  Potentials 
Using  the  Old  and  New  Schemes 
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DSX 


Old  Script 


?  C 


r8 


Charge_Suifaces 
A  xmlns 
A  Prefix 
c  ReadObject 
C  OpenOatabase 
C  SetConductorBias 
C  ComputeTransverseCurrents 
C  SetimtialConductorPotentiai 
c  SetlmtialPotenbais 
C  WntePotenbals 
Embed_ObjectJn_Grid 
PotentialsJn_Space 
Create_Pai  tides 
Loop 

A  Iterations 
A  StartAl 
c  Track_Partides 

C  Charge_Surfaces 

A  xmlns 

A  Prefix 

o-  C  OpenDatabase 
C  InitializeCalculabon 
o-  C  SetHlummabon 
o-  C  SetEnvironment 
O-  C  SetCondudorBias 

C  ComputeTransverseCurrents 
C  UseTrackedCurrents 
C  ZeroCurDerivAlgontfim 
c  PrepareChargeMatnx 
o-  C  DoTrackTimeStep 
c  Potentials_in_Space 
c  Create_Partides 
C  Save_Files 


New  Script 


Script  - 


-  C  Read_Objed 

-  C  Emoed_Objed_in_Gnd 

-  C  Potential  s_in_Space 

-  C  Create_Parbdes 

-  C  Track_Paitides 

9  C  Charge_Surfaces 
»•  O  Settllummation 
o-  C  SetEnvironment 
o-  C  SedmtialCondudorPotential 
*-  O  SetCondudorBias 
C  UseTrackedCurrents 
®-  Q  SetParameters 
C  DoTrackTimeStep 
9  C  Loop 

A  Iterations 
A  StartAl 

*•  C  Potentiais_in_Space 
®-  C  Create_Paitides 
O-  C  Track_Parddes 
9  C  Charge_Suifaces 
o-  Q  Settllummation 
>■  C  SetEnvironment 
C  SetCondudorBias 
c  UseTrackedCurrents 
O-  C  SetParameters 
C  DoTrackTimeStep 
o-  C  Save_Flles 


s3 


30 

1 


Figure  4.  Comparison  of  Scripts  Used  to  Specify  the  DSX  Charging  Problem  Using  the  Old  and  New  Schemes 

Nascap-2k  standard  test  problems  were  used  to  verify  that  the  scripting  changes  are  adequate  and 
that  they  were  correctly  implemented.  The  suite  of  standard  test  problems  was  expanded  to 
include  at  least  one  computation  that  uses  each  of  the  available  script  commands.  Using  the 
existing  scripts,  the  standard  test  problems  were  executed  using  the  code  from  before  the  changes 
and  using  the  code  from  after  the  changes.  All  the  differences  in  the  text  output  files  were 
examined  and  the  causes  determined.  All  issues  were  fixed. 

For  the  test  problems  that  use  a  default  script,  a  new  default  script  was  generated.  For  test 
problems  that  use  a  modified  script,  the  script  was  recreated  using  the  new  formulation.  The 
calculations  were  done  using  the  new  script  and  all  differences  in  text  output  examined  and  the 
causes  determined.  All  issues  were  fixed. 

3.  GENERALIZATION  OF  NASCAP-2K  CHARGING  CAPABILITIES 

Nascap-2k  was  modified  to  allow  for  a  wider  variety  of  charging  calculations  to  be  easily 
specified.  Table  1  lists  the  charging  calculations  that  can  now  be  specified  within  the  user 
interface. 
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Table  1.  Charging  Calculations  Allowed  by  Default  (Blue  Text  Indicates  New  Options) 


Geosynchronous 

•  Analytic  currents  (Maxwellian,  Double  Maxwellian,  Kappa,  Tabular) 

•  No  velocity,  Only  H+  ions 

LEO  or  Plume 

•  Analytic  currents  (Convected  Maxwellian  for  user  specified  ion  species  and  for  electrons) 

•  Tracked  particle  currents 

•  Tracked  ion  and  analytic  electron  currents  (Convected  Maxwellian) 

Auroral 

•  Analytic  currents  (Convected  Maxwellian  for  user  specified  ions  species  and  Fontheim  for 
electrons) 

•  Tracked  ion  and  analytic  electron  currents  (Fontheim) 

Interplanetary 

•  Analytic  currents  (Convected  Maxwellian  for  user  specified  ions  species  and  for  electrons) 

•  Tracked  ion  and  analytic  electron  currents  (Convected  Maxwellian) 

LEO/Time  dependent  plasma 

•  Tracked  particle  currents 

•  Tracked  ion  and  analytic  electron  currents  (Convected  Maxwellian) 


The  user  interface  was  modified  in  two  ways.  The  code  that  enables  and  disables  the  options 
available  on  the  Problem  tab  was  modified  to  allow  for  the  wider  variety  of  calculations  to  be 
specified.  Also,  the  code  that  writes  the  default  scripts  was  modified  to  allow  for  the  additional 
calculations.  At  the  same  time  the  default  script  code  was  streamlined. 

The  most  important  changes  made  are  those  in  the  code  that  computes  the  current  to  a  surface 
and  its  derivative.  The  code  was  both  modified  and  cleaned  up.  The  tight  coupling  between  the 
type  of  environment  and  the  current  calculation  was  reduced.  The  code  that  computes  the  current 
and  the  code  that  computes  its  derivative  were  reorganized  to  make  their  structure  clearer  and 
extensive  comments  were  added. 

Analytic  formulas  for  the  derivative  of  the  current  with  respect  to  the  potential  for  the  Fontheim 
and  Convected  Maxwellian  distributions  were  computed.  For  the  Fontheim  distribution  the 
formula  is  exact,  and  for  the  Convected  Maxwellian  it  is  approximate.  These  formulas  were 
implemented  in  the  code  to  replace  the  numeric  derivatives  previously  used. 

The  capability  to  compute  and  use  secondary  electron  currents  in  all  charging  calculations  with 
tracked  currents  was  added.  The  yield  is  computed  using  analytic  currents  with  the  correct  ion 
masses. 

The  current  to  surfaces  was  computed  in  four  separate  methods  (subroutines).  The  code  was 
reorganized  so  that  all  four  methods  call  a  single  method  containing  the  current  computation 
code. 
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Sample  calculations  for  charging  (current  balance)  in  a  cold,  dense  plasma  indicated  that  the  last 
step  in  the  computation  of  the  derivative  of  the  current  with  respect  to  the  potential  needed  to  be 
changed.  Previously,  the  derivative  was  always  set  to  be  large  enough  so  that  the  estimated 
change  in  potential  during  the  next  timestep  is  at  least  1  V.  For  stability,  the  minimum  estimated 
change  in  potential  needed  to  be  reduced. 

The  following  example  was  used  to  test  the  new  capabilities:  charging  in  a  cold,  dense  plasma 
using  currents  computed  for  a  convected  Maxwellian  plasma.  The  Object  Toolkit  model  was 
constructed  at  AFRL  to  model  a  LEO  spacecraft.  The  geometry  is  shown  in  Figure  5,  with  the 
results  in  Figure  6.  The  solar  array  is  divided  into  four  sections  with  different  potentials  from 
chassis  plus  50  V  to  chassis  plus  200  V.  The  potentials  on  the  gold  surfaces  reflect  these  values. 
The  results  are  as  expected.  The  chassis  potential  rapidly  decreases  to  -140  V.  The  insulating 
OSR  and  solar  cell  surfaces  show  immediate  potential  swings  with  the  underlying  conductors. 
The  insulating  surfaces  then  differentially  charge  to  small  potentials  reflective  of  the  plasma 
temperature. 


Materials 


Graphite 


Gold 


Figure  5.  Geometry  and  Materials  for  First  Test  Case  for  Charging  in  LEO  with  Analytic  Currents  from  a 

Convected  Maxwellian  Distribution 
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Figure  6.  Potentials  Resulting  from  Charging  Calculation  using  Geometry  of  Figure  5 

4.  COMPREHENSIVE  TESTING 

Comprehensive  testing  of  the  new  code  was  performed. 

Calculations  done  earlier  (both  the  standard  test  problems  and  others)  were  repeated  with  the 
new  code  and  either  the  old  or  regenerated  scripts.  To  speed  up  these  comparisons,  flags  were 
modified  to  reduce  non-meaningful  differences  in  text  files  (such  as  the  time  since  execution 
started).  Also,  Nascap-2k  was  modified  so  that  random  number  generation  can  occur  with  a  fixed 
seed  or  with  a  random  seed.  Test  problems  that  use  the  random  number  generator  were  adjusted 
to  use  a  fixed  seed. 

Also,  calculations  enabled  by  recent  code  enhancements  were  performed  and  the  results 
compared  with  expected  values. 

Also,  the  standard  suite  of  test  problems  was  executed  on  LINUX  and  the  appropriate 
comparisons  done. 

Examples  of  this  testing  follow. 

4.1.  Big  Grid  Test 

A  “Potentials  in  Space”  calculation  was  done  with  a  grid  with  a  large  number  of  volume 
elements.  The  computational  grid  is  121  x  101x51  with  mesh  spacing  of  6.7  cm.  The  potentials 
about  a  -10  V  parallelepiped  in  a  4  cm  Debye  length  plasma  were  computed.  The  results  are 
shown  in  Figure  7. 
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4.2.  Magnetic  Dipole 


Verified  that  the  capability  to  track  particles  in  a  dipolar  magnetic  field  still  functions  correctly 
by  tracking  electrons  from  the  edge  of  a  sheath  to  a  +10  V,  1  m  radius  sphere  with  a  300  Am2 
magnetic  dipole  at  its  center.  The  potentials  and  the  trajectories  are  shown  in  Figure  8. 


Figure  8.  Electron  Trajectories  toward  a  +10  V  Sphere  with  a  Magnetic  Dipole  at  Its  Center  (Potentials  in  a 
Plane  through  the  Center  of  the  Sphere  Are  Also  Shown) 

4.3.  Simple  vxB  Test 

An  earlier  calculation  done  while  developing  the  capability  to  include  v  x  B  potentials  in 
charging  calculations  was  repeated  to  verify  that  these  potentials  are  still  included  correctly  after 
the  recent  changes  of  the  scripting  and  charging  coding.  The  existing  (previously  developed) 
script  (shown  in  Figure  9)  was  used.  The  geometry  used  is  a  metal  cylinder  with  Teflon  collars  as 
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shown  in  Figure  10.  The  magnetic  field  is  taken  to  be  7.1  gauss,  and  the  velocity  is  taken  to  be 
600  km/s,  for  an  induced  potential  of  1700  V  along  the  4  m  long  cylinder.  The  plasma  density 
starts  at  10  m'  and  increases  by  one  order  of  magnitude  each  second  to  10  m'  from  4  seconds 
to  the  end  of  the  calculation  at  10  seconds.  The  electron  and  ion  temperatures  are  both  taken  as  1 
eV.  A  constant  timestep  of  0.2  seconds  was  used. 

The  results  shown  in  Figure  1 1  through  Figure  13  are  as  expected.  The  v  x  B  potentials  appear 
on  the  exposed  conductor  surface  elements.  The  insulators  are  initially  at  the  same  potential  as 
the  underlying  conductor  (zero  differential  potential)  and  charge  to  zero  absolute  potential,  at 
first  slowly  and  then  more  rapidly  as  the  plasma  density  is  increased,  as  shown  in  Figure  12.  The 
charging  current  to  insulators  at  the  negative  end  shows  roughly  an  order  of  magnitude  increase 
each  time  the  plasma  density  is  increased,  while  the  charging  current  to  insulators  at  the  positive 
end  remains  small,  as  shown  in  Figure  13. 


<?xml  version  ="1.0"  encoding="UTF-8n  standalone=nno"  ?> 

-  <SCRIPT  > 

-  cCOMMAND  Prefix="vxBDiel"  cmd="Charge_Surfaces"> 

<COMMAND  FileName="vxBDielObject.xml"  cmd="ReadObject"  /> 

<COMMAND  Prefix  ="vxBDiel"  cmd-  OpenDatabase  /> 

<COMMAND  cmd="InitializeCalculation"  > 

cCOMMAND  Value-'O.O"  cmd="SetIllumination"  x-'l.O"  y-’O.O"  z="0.0"  /> 

-  cCOMMAND  cmd="SetEnvironment":> 

cEnvironment  begintime-'O.O"  nel-'lOOOOOO.O"  nil="1000000.0"  tel-'l.O"  til="1.0"  type-’GEO"  /> 
cEnvironment  begintime-’l.O"  fraction l=n1.0''  massl="2.6568624E-26"  nel="1.0E7"  nil="1.0E7"  numSpecies="l" 
shadowlons-'yes"  tel-'l.O"  til-'l.O"  type-'GEO"  /> 

cEnvironment  begintime="2.0"  fraction  1  ="1.0"  massl="2.6568624E-26"  nel="1.0E8"  nil="1.0E8"  numSpecies="l" 
shadowlons-’yes"  tel-'l.O"  til-'l.O"  type-'GEO"  /> 

cEnvironment  begintime="3.0"  fractionl-'l.O"  massl="2.6568624E-26"  nel="1.0E9"  nil="1.0E9"  numSpecies="l" 
shadowlons-'yes"  tel-'l.O"  til="1.0"  type-'GEO"  /> 

cEnvironment  begintime="4.0"  fractionl-'l.O"  massl="2.6568624E-26"  nel="1.0E10"  nil="1.0E10"  numSpecies="l" 
shadowIons="yes"  tel="1.0"  til="1.0"  type="GEO"  /> 
c/COMMAND  > 

cCOMMAND  cmd=  SetVelocity  x=  600000.0  y="0.0"  z='0.0  /> 
cCOMMAND  cmd="SetBField"  x="0.0"  y="7.1E-4"  z="0.0"  /> 
cCOMMAND  Index  ="1"  Value-'O.O"  cmd="SetInitialConductorPotential"  /> 
cCOMMAND  cmd="SetInitialPotentialsM  /> 
cCOMMAND  Value  ="0.0"  cmd=  SetVXBPotentials  > 
cCOMMAND  cmd="PrepareChargeMatrix"  /> 

-  cCOMMAND  CloseFiles="True"  cmd="DoTimeSteps"> 

cTimeParams  begintime-'O.O"  endtime-'lO.O"  maxdt="0.2"  mindt="0.2"  nsteps="50"  /> 

^COMMAND  > 
c/COMMAND  > 

cCOMMAND  InputFileName  ="vxBDiel_n2kdyn_in.txt"  OutputFileName  ="vxBDiel_n 2kdyn_out.txt" 
cmd="Embed_Object_in_Grid"  /> 

cCOMMAND  InputFileName -'vxBDiel_potent_0_in.txt"  Iteration="0"  0utputFileNarne-  vxBDiel_potent_0_out.txt" 
cmd="Potentials_in_Space"  /> 
c/SCRIPT  > 


Figure  9.  Script  Used  for  “Simple  yxB”  Test 
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Figure  10.  Object  Toolkit  Object  Used  in  Verification  Calculation 


Figure  11.  Final  (Steady  State)  Surface  Potentials  for  the  Verification  Calculation 
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Potential 


Differential  Potential 


Figure  12.  Absolute  and  Differential  Potential  as  a  Function  of  Time  (and  Density)  at  Selected  Insulating 
Surface  Elements  (Two  Surfaces  from  Each  Insulating  Collar  are  Shown) 

Charging  Current 


Figure  13.  Charging  Current  to  Insulating  Elements  at  Each  End  of  the  Cylinder  as  a  Function  of  Time  (and 
Density) 

4.4.  Juno 

An  earlier  calculation  of  the  potentials  expected  on  the  Juno  spacecraft  during  a  polar  pass  over 
the  Jupiter  pole  for  specified  conditions  was  repeated  to  verify  that  the  capabilities  used  for  that 
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calculation  still  work  correctly.  The  calculation  tests  the  ability  to  simulate  a  spinning  satellite 
moving  rapidly  through  a  high  magnetic  field,  time-dependent  environment.  In  particular,  we 
wanted  to  be  sure  that  the  differential  potentials  are  preserved  across  successive  Charge 
Surfaces  commands,  and  that  the  time-dependent  environment  was  updating  properly, 

The  previously  generated  script  was  used,  with  the  exception  that  “starttime”  keyword  for  the 
environment  were  changed  to  “begintime”  keywords  to  match  the  new  convention.  The 
spacecraft  rotates  at  5  rpm.  The  velocity  vector  (60  km/s)  lies  in  the  plane  of  rotation,  along  with 
a  1  gauss  magnetic  field  perpendicular  to  the  velocity.  The  magnetic  field  normal  to  the  plane  of 
rotation  was  10  gauss,  giving  a  v><B  field  in  the  plane  of  rotation  of  60  Y/m.  Surface  potentials 
on  the  bottom  (conductive)  side  and  on  the  top  (solar  cell)  side  are  shown  in  Figure  14. 

The  code  correctly  executed  the  test  case.  It  was  noted  that  the  first  environment  in  the  first  set  is 
automatically  reset  to  have  a  “begintime”  of  0.0. 


Figure  14.  Surface  Potentials  on  Juno  for  the  Bottom  (Conductive)  Side  (Left)  and  the  Top  (Insulating  Solar 
Cell)  Side  (Right). 

4.5.  SPEAR  3 

Calculation  of  the  potentials  about  the  SPEAR  III  rocket  and  its  magnetically  limited  current 
collection  was  performed.  This  calculation  had  not  been  done  since  about  1995.  The  Patran 
model  of  SPEAR  III  was  imported  into  Object  Toolkit  and  materials  and  conductor  numbers 
assigned  to  surfaces.  The  density  was  set  to  1011  m'3,  the  temperature  to  0.3  eV,  and  the  ion 
species  to  Oxygen.  Bias  of  10  kV  was  applied  to  the  sphere,  with  750  V  increments  on  the 
graded  boom,  and  the  body  was  set  to  -3  kV.  Potentials  were  calculated  using  the  “Non-linear” 
formulation,  giving  the  results  shown  in  Figure  15.  Then  particles  were  tracked  with  no  magnetic 
field  and  with  0.3  gauss  along  the  boom  or  normal  to  the  plane  of  the  model.  The  resulting 
currents,  positive  to  the  body  and  negative  to  the  boom  and  sphere,  are  shown  in  Table  2.  They 
follow  the  expected  pattern. 
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Figure  15.  Potentials  about  the  SPEAR  III  Model,  with  Parameters  as  Given  in  Text 


Table  2.  Calculated  Currents  to  the  SPEAR  III  Model  for  Parameters  Given  in  Text  and  for  Three  Magnetic 

Field  Configurations 


Chassis 

potential 

Sphere 

potentials 

Magnetic  field 
(gauss) 

Current  to 
body  (A) 

Current  to 
sphere (A) 

Current  to 
boom  (A) 

Total 

current (A) 

-3000 

7000 

(0,0,0) 

0.00103 

-0.0306 

0 

-0.030 

-3000 

7000 

(0, 0.3,0) 

0.00103 

-0.0108 

-0.00032 

-0.010 

-3000 

7000 

(0,0, 0.3) 

0.00103 

-0.0085 

-0.00837 

-0.016 

4.6.  MMS 

Three  earlier  calculations  of  the  charging  of  the  MMS  spacecraft  were  repeated  to  verify  that  a 
standard  long  Debye  length  charging  calculation  is  still  easily  defined  and  that  the  results  remain 
as  earlier.  The  three  calculations  have  different  environments  and  timesteps.  All  three  were 
executed  with  the  old  script  and  with  a  rebuilt  default  script.  The  results  are  identical  within 
expected  numeric  variation. 

4.7.  Solar  Probe  Plus 

An  earlier  calculation  of  the  charging  of  the  Solar  Probe  spacecraft  was  repeated  to  verify  that 
the  functionality  used  to  compute  charging  in  the  Solar  Wind  continues  to  work  correctly.  Once 
the  parameters  were  adjusted  to  be  identical  to  those  used  previously,  the  same  results  were 
obtained,  within  the  expected  numeric  variation. 
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4.8.  Emitter  Test  Case 


In  order  to  explore  and  test  Netscape's  ability  to  model  an  emitter,  the  floating  potential  of  a 
cylinder  with  an  emitter  on  its  side  was  computed.  The  current  balance  is  between  the  emitter 
electron  current  and  sheath  electrons  collected  in  Parker-Murphy  fashion  from  a  LEO  plasma. 

A  simplified  version  of  this  test  case  has  been  added  to  the  Nascap-2k  User  s  Manual  as  an 
example. 

4.8.1.  The  Parameters 

The  object  is  a  1  m  long,  0.4  m  diameter,  aluminum  cylinder,  as  shown  in  Figure  16.  A  current  of 
0.3  Am'2  of  electrons  is  emitted  from  a  surface  with  area  0.0345 1  m'2,  for  a  total  current  of 
10.35  mA.  The  environment  parameters  are  set  as  shown  in  Figure  17.  The  magnetic  field  is 
parallel  to  the  cylinder  axis,  and  an  extra  species  is  defined  for  emitted  electrons.  On  the 
Problem  tab  the  environment  is  specified  as  “LEO  or  Plume”  and  the  problem  type  is  specified 
as  “Charging”  with  “Tracked  Particle  Currents”.  The  conductor  is  specified  as  floating  with  an 
initial  potential  of  +5  V.  The  computational  space  is  three  nested  grids  with  a  20  cm  resolution 
outer  grid  and  a  5  cm  resolution  inner  grid.  The  emitter  properties  are  shown  in  Figure  18. 


Figure  16.  Aluminum  Cylinder  Showing  Emitter  Surface  Element 
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Figure  17.  Environment  Parameters 


I  Emitter  EGun  Properties 


CurrentDensity 

NumLocations 

NumPhis 

NumThetas 

MinTheta 

MaxTheta 

NumEnergies 

MinEnergy 

MaxEnergy 


2999.0 


3001.0 


OK 


Figure  18.  Emitter  Parameters 


4.8.2.  Building  the  Script 

The  script  for  this  problem  is  nonstandard.  The  following  procedure  was  used. 

1 .  On  the  Charging  tab,  set  20  0. 1  ms  timesteps. 

2.  On  the  Particles  tab,  make  the  appropriate  selections  for  the  emitter  electrons,  as  shown 
in  Figure  19. 

3.  On  the  Script  tab,  build  the  script.  Make  sure  the  “Automatically  overwrite  files”  box  is 
unchecked. 
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4.  For  all  the  Create  Particles  and  Track  Particles  commands,  set  the  index  of  the  input 
file  to  0,  as  shown  in  Figure  20. 

5.  Click  “Save  Files.”  Click  “Yes  to  All”  on  the  overwrite  question. 

6.  On  the  Particles  tab,  make  the  appropriate  selections  for  the  sheath  electrons,  as  shown 
in  Figure  21. 

7.  Return  to  the  Script  tab  and  rebuild  the  script. 

8.  Make  the  same  edit  specified  in  item  4  above. 

9.  Save  the  files  with  overwrite,  as  in  item  5  above. 

10.  Duplicate  each  of  the  three  Create_Particles  commands,  as  shown  in  Figure  22. 

1 1.  In  one  of  each  pair  of  Create_Particles  commands,  change  “Electron”  to 
“ElectronsEmitted”  (three  places)  as  shown  in  Figure  23.  (The  Track_mode, 
Creation_mode,  and  Species  values  are  used  to  write  the  input  file.  When  existing  input 
files  are  used  (as  in  step  13  below)  the  values  are  ignored.) 

12.  Edit  the  Tracker  input  file  to  specify  both  species,  as  shown  in  Figure  24. 

13.  Run  the  script,  this  time  answering  “No  to  All”  on  the  overwrite  question. 

Note  that  several  methods  are  available  to  construct  the  desired  input  files  and  script. 

Charged  particles  used  for  surface  current  calculation 

■  Initial  Particle  Distribution - 

O  Sheath  Potential  Value  (V): 

C  B  Field  Number  per  zone: 

O  Boundary  * 

O  Charge  Exchange  Cross  Section(10'20  Mb0-0  1 

<•)  Surface  *  |  Emitter  |  ▼]  EGun  ▼j  Current  Density  (Am'2)j0.300  | 

O  None  of  above 

G  External  File  Filename:  |  Browse  | 

’Additional  parameters  on  advanced  screen. 

■  Particle  Species - 

■  Electron 

I  Oxygen _ 

■  ElectronsEmitted 


Advanced 


Figure  19.  Particles  Tab  Specification  for  Emitter  Electrons 
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Figure  20.  Editing  of  Script  to  Simplify  Electron  Emitter  Files 
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Charged  particles  used  for  surface  current  calculation 


■  Initial  Particle  Distribution - 

®  Sheath 

Potential  Value  (V):|o.300 

O  B  Field 

Number  per  zoneqO 

O  Boundary  * 

O  Charge  Exchange 

O  Surface* 

Current  Density  (Am‘2)j0.300 

O  None  of  above 

□  External  File  Filename: 

Browse 

"Additional  parameters  on  advanced  screen. 


-  Particle  Species  — 

■  Electron 

■  Oxygen 

■  Electron  sEmitted 


Advanced 


Figure  21.  Particles  Tab  Specification  for  Sheath  Electrons 


Figure  22.  Script  with  Duplicated  Create  Particles  Commands 
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Figure  23.  Editing  of  a  Pair  of  Create  Particles  Commands  to  Create  both  Emitter  and  Sheath  Electrons 


4  Emitter_tracker_trajE_0_in.txt  -  Notepad 

File  Edit  Format  View  Help 

Comment  —  Tracker  Input  File  for  Mode  Surface  Currents 
Commen  t 

Comment  —  File  Prefix  ... 

PREFIX  Em  it  ter 


Comment 

Comment  --  Options 

PROCESS 

UP  DATE  TIME 

CUT_PLANE 

TRACK_  TIME 

DX__MAX 

MAX-STEP 

SA  VE _INTER  VAL 

B_FIELD 

MIX 

SPLIT 

SPECIES 

SPECIES 

TRACK-EL  EC  TRONS 


TRAJECTORIES 

FALSE 

X  0. 0000E00 

1 . 0000E 00  seconds 

1 . 0000E-01  loca 7  gr i d  un i ts 

9999 

1  STARTING  1 

0.0000E00  0.  0000 E 00  3. 0000E-05 
0. 0000E00 
Off 

E lectron 

E lectronsEmi tted 

Save  charge  density  to  electron 


arrays 

Comment 


Comment  --  miscellaneous 


X_ LIMIT 

Y_ LIMIT 

Z— LIMIT 

X—PL  OT—L IMIT 

Y—PLOT—L IMIT 

Z—PL  OT—L  IMIT 

DIAG-TRACKER 

DIAG_DYNALIB 

TIMER 

END 


0. 0000E00  0. 0000E00 
0.0000E00  0.  0000E00 
0. 0000E00  0.  0000E00 
0.0000E00  0.  0000E00 
0. 0000E00  0. 0000 E 00 
0.0000E00  0.  0000 E 00 
1 
1 
1 


Figure  24.  Edit  of  Tracker  Input  file.  Note  That  Both  Species  (Sheath  Electrons  and  Emitter  Electrons)  are 

Specified  for  Tracking 


4.8.3.  Results 

The  resulting  potential  versus  time  is  shown  in  Figure  25.  Early  on,  the  charging  rate  is  limited 
by  the  dJ/dV,  which  is  set  to  0.75J/V.  The  potential  increases  until  the  sheath  current  balances  the 
emitter  current  at  about  840  V.  The  charging  rate  continues  to  be  limited  because  the  code  sets 
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.  i.e,  dl/dV  is  large  and  negative  even  though  I  is  near 


dl 

~dV 


X 

elem 


dJ_ 

~dV 


-0.75X 


elem 


J 

«  -0.75 

/ 

V 

V 

zero  due  to  cancellation  of  the  emitted  and  sheath  electrons. 


The  net  current  versus  potential  is  shown  in  Figure  26.  At  about  840  V  the  sheath  electron 
current  cancels  the  emitter  current.  Note  that  the  net  current  is  obtained  by  taking  the  average 
“Tracked  Electron  Current”  to  all  surface  elements  and  multiplying  by  the  total  area  of  1.4823 
m2. 


If  the  option  is  used  to  set  the  stabilizing  current  derivatives  to  zero,  then  the  correct  (using  the 
vacuum  capacitance)  charging  timescale  is  achieved,  as  shown  in  Figure  27. 

Potential 


Figure  25.  Potential  versus  Time  for  the  Emitter  Problem 


Figure  26.  Net  Current  versus  Potential. 
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Figure  27.  Charging  Dynamics  with  the  Stabilizing  Current  Derivatives  Set  to  Zero 
4.9.  Detector  Test  Case 

In  order  to  explore  and  test  Nascap-2k’s  ability  to  model  a  detector,  the  current  to  a  detector  with 
a  guard  ring  was  computed. 

This  test  case  has  been  added  to  the  Nascap-2k  User  s  Manual  as  an  example. 

4.9.1.  Object 

The  object,  intended  to  represent  a  detector  with  a  rudimentary  collimator,  is  shown  in  Figure  28. 
The  detector  (cyan  surface  element)  and  the  surrounding  eight  surface  elements  (guard  ring)  are 
defined  as  conductor  3,  and  are  biased  +5  V  relative  to  the  remainder  of  the  plate  and  the 
collimator  (yellow  surface  elements). The  definition  of  the  detector  is  shown  in  Figure  29.  It 
specifies  the  emission  of  test  particles  from  four  locations  on  the  detector  in  144  directions  and 
with  20  energies  ranging  from  5  eV  (0  eV  of  total  energy)  to  7  eV  (2  eV  of  total  energy). 

4.9.2.  Problem  Tab 

On  the  Problem  tab,  the  environment  is  specified  as  “LEO  or  Plume”  and  the  problem  type  as 

“Potentials  in  Space  or  Detector  Analysis”  with  “Analytic  Space  Charge”  and  “Surface 
Currents”. 
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Figure  28.  Detector  Test  Case  Object. 


i  -.i  Detector  PartDetect  Properties  .  1  1  ©“J  l  «=  1  0 

NumLocations 

NumPhis 

NumThetas 

MinTheta 

MaxTheta 

NumEnergies 

MinEnergy 

MaxEnergy 

4 

12 

12 

0.0 

1.57 

20 

5.0 

7.0 

1 

OK  |  |  Cancel 

Figure  29.  Definition  of  the  Detector  Properties 


4.9.3.  Environment  Tab 

11  3 

The  Environment  tab  specifies  a  plasma  with  density  10  m'  ,  temperature  0.3  eV,  and  no 
motion  or  magnetic  field. 

4.9.4.  Applied  Potentials  Tab 

On  the  Applied  Potentials  tab,  conductor  1  is  fixed  at  0  V  (ground),  and  conductor  2 
(collimator)  and  conductor  3  (detector  and  guard  ring)  are  biased  at  0  and  +5  V  respectively. 

4.9.5.  Space  Potentials  Tab 

On  the  Space  Potentials  tab,  the  “Non-linear”  potential  formulation  is  chosen. 
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4.9.6.  Particles  Tab 


The  Particles  tab  (Figure  30)  specifies  electron  emission  from  a  surface  using  the  “Detector” 
treatment,  with  specifications  for  the  “PartDetect”  detector  that  was  defined  during  object 
definition  by  Object  Toolkit. 


Charged  particles  used  for  surface  current  calculation 


initial  Particle  Distribution - 

G  Sheath 
O  B  Field 
G  Boundary  * 

O  Charge  Exchange 
®  Surface*  | Detector  | 

O  None  of  above 
□  External  File  Filename: 


PartDetect 


Potential  Value  PQql.OOO 
Number  per  zone: 


‘Additional  parameters  on  advanced  screen. 

-  Particle  Species - 

■  Electron 


Advanced 


Figure  30.  Particles  Tab  for  the  Detector  Simulation 


4.9.7.  Script  Tab 

With  the  above  settings,  the  default  script,  Figure  31,  is  used.  Run  the  script  to  obtain  the  results 
discussed  in  the  following  sections. 


Run  Script 


Edit  Script 


o  C  Read_Object 

C  Embed_Object_in_Grid 
o  C  Potentials_in_Space 
o-  G  Create_Particles 
o  C  Track_P articles 


Figure  31.  Script  for  the  Detector  Problem  (Default  Script) 


4.9.8.  Results  Tab 

With  no  charging  and  only  one  surface  element  with  current,  there  is  very  little  information  to  be 
viewed  on  the  Results  tab,  seen  in  Figure  32.  It  is  possible  to  see  that  one  surface  element  (the 
detector  surface  element)  is  assigned  a  current  density  of  -3.257  mAm"2. 
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Figure  32.  Results  Tab  Following  Execution  of  the  Detector  Problem 

4.9.9.  Results  3D  Tab 

The  Results  3D  tab,  shown  in  Figure  33  is  used  to  obtain  the  results.  The  properties  of  the 
detector  surface  element  confirm  the  current  density  noted  above.  Multiplying  by  the  surface 
element  area  gives  4.0  pA  as  the  current  to  the  detector.  The  potential  bows  out  significantly 
from  the  collimator  aperture,  suggesting  that  the  current  is  somewhat  greater  than  the  aperture 
area  times  the  electron  thermal  current.  Since  the  aperture  approximates  a  circle  of  2.5  cm  radius 
and  the  electron  thermal  current  is  1.468  mA  nT  ,  the  planar  current  through  the  aperture  would 
be  only  2.9  pA.  (By  comparison,  a  sheath  calculation  with  the  sheath  potential  set  to  0.2  V  gives 
3.4  pAto  the  detector  surface  element.) 

4.9.10.  Tracker  Output 

The  Tracker  output,  the  relevant  excerpt  of  which  is  shown  in  Figure  34,  sheds  additional  light 
on  the  calculation.  Of  11520  test  particles  launched  from  the  detector  surface  element,  1365  left 
the  primary  grid  (meaning  that  they  represent  electrons  from  the  environment),  while  10155 
struck  object  surfaces  (meaning  they  represent  phase  space  that  does  not  connect  to  the 
environment).  Applying  environmental  factors  to  the  “lost”  current  particle  weights  yields  the  4.3 
pA  noted  above,  and  is  assigned  to  conductor  3  “alum”  surface  elements  (representing  the 
detector).  Applying  environmental  factors  to  the  10155  test  particles  that  struck  the  object  yields 
5.14  pA;  this  is  not  necessarily  a  meaningful  number  as  these  particles  do  not  represent  electrons 
from  the  environment,  but  gives  a  better  estimate  than  the  raw  particle  numbers  or  weights  of 
how  much  of  the  detector’s  phase  space  is  blocked. 
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Figure  33.  Results  3D  Tab,  Showing  Properties  of  Detector  Surface  Element,  and  Showing  Bowing  Out  of 

Potentials  through  the  Collimator 


|  Detector_tracker_trajE_0_out.txt  -  Notepad 

looi^ji^iiHis^r 

File  Edit  Format  View  Help 

Trackr:  total  of  11520  new  particles . 

0  were  partially  tracked. 
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1365  went  off  primary  grid. 
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Proces:  current  to  object:  -5 . 1393E-06  amps. 

lost  current (of f  grid):  -4 .  3065E-06  amps, 
trapped  current  :  0.  OOOOE+OO  amps, 

other  current  :  0.  OOOOE+OO  amps. 


Cond. 

kapt 

tefl 

a  lum 

gold 

osr 

b  lac 

so  la 

grap 

Tota  1 

1 

0.  OE+OO 

0.  OE+OO 

0.  OE+OO 

0.  OE+OO 

0.  OE+OO 

0.  OE+OO 

0.  OE+OO 

0.  OE+OO 

0.  OOOE+OO 

2 

0.  OE+OO 

0.  OE+OO 

0.  OE+OO 

0.  OE+OO 

0.  OE+OO 

0.  OE+OO 

0.  OE+OO 

0.  OE+OO 

0.  OOOE+OO 

3 

0.  OE+OO 

0.  OE+OO 

-4. 3E-06 

0.  OE+OO 

0.  OE+OO 

0.  OE+OO 

0.  OE+OO 

0.  OE+OO 

-4. 307 E -06 

Tot  a  1 

0.  OE+OO 

0.  OE+OO 

-4. 3E-06 

0.  OE+OO 

0.  OE+OO 

0. OE+OO 

0.  OE+OO 

0.  OE+OO 

-4.  307 E -06 

Figure  34.  Tracker  Output  from  Detector  Run 
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5.  MEMO  SUMMARIZING  RESOURCES  FOR  INTERNAL  CHARGING 
COMPUTATIONS 


leidos 


Date:  January  23,  2014 

To:  Dale  Ferguson  and  Adrian  Wheelock 

From:  Victoria  Davis  and  Myron  Mandell 

Subject:  Summary  of  Resources  Regarding  Internal  Charging  Environment  and  Internal 
Charging  Computation 

The  following  is  a  compilation  of  resources  we  have  used  at  times  in  assessing  risk  due  to 
Internal  Charging. 

Proposed  Worst  Case  and  Design  Environments 
Fennell  Environment 

This  environment  is  described  in  a  paper  by  Joe  Fennell  [Fennell  et  al.,  2000].  The  paper 
proposes  a  worst  case  spectrum  for  geosynchronous  orbit  based  on  measurements  taken  after  the 
March  24,  1991  magnetic  storm  event.  They  also  propose  a  HEO/Molniya  worst  case  spectrum. 

Flumic 

Flumic  was  developed  at  QinetiQ  and  is  documented  in  two  papers,  [Wrenn,  et  al.,  2000]  and 
[Rodgers  et  al.,  2003].  The  model  has  been  validated  by  a  comparison  with  flight  data  [Evans,  et 
al,  2007]  and  [Taylor,  et  al,  2009].  Flumic  is  included  in  Spenvis  for  use  with  DICTAT.  Victoria 
has  the  formulas  implemented  in  an  Excel  spreadsheet.  The  Flumic  spectrum  depends  on  L,  solar 
cycle,  and  season.  The  inner  belt  model  also  has  a  B/Bo  dependence.  In  general,  it  gives  more 
internal  charging  than  other  models  due  to  the  spectral  shape. 

AES 

AE8  is  a  climatological  average  of  the  environment.  As  such,  it  is  useful  for  estimating 
degradation  of  solar  cells  over  a  period  of  years.  It  is  sometimes  used  to  estimate  internal 
charging  after  enhancement  by  a  factor  on  the  order  of  ten. 

NASA-HDBK-4002 

The  Handbook  has  a  figure  showing  a  proposed  worst  case  environment  for  internal  charging  on 
geosynchronous  spacecraft.  At  one  time,  Victoria  read  the  numbers  off  the  figure  and  has  done 
calculations  using  it. 

Frederickson  fit 

Robb  Frederickson  had  the  differential  fluence  over  the  CRRES  five-hour  half-orbits  fit  to  a 
functional  form.  These  fits  are  tabulated  in  a  NASA  report  [Frederickson  and  Brautigam,  2004]. 
Victoria  determined  that  the  orbit  600  spectrum  gives  the  most  internal  charging. 
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Software  to  Compute  Deep  Dielectric  Charging 
NUMIT 

Several  versions  of  Robb  Frederickson’ s  Fortran  code  NUMIT  exist.  Joe  Minow,  Insoo  Jun,  and 
Brian  Beecken  have  all  evaluated  the  code  and  made  modifications.  I  believe  that  David  Cooke 
has  Brian  Beecken  on  contract  to  create  a  standard  version.  Insoo  Jun  announced  a  similar  effort 
recently. 

SAIC  Deep  Dielectric  Charging  Code 

The  Leidos  proprietary  code  is  a  more  flexible  stand-alone  version  of  the  code  originally 
implemented  in  the  SEE  Spacecraft  Charging  Flandbook.  It  implements  the  same  algorithms  as 
NUMIT  and  Robb  was  heavily  involved  in  the  original  implementation.  It  has  been  verified  to 
give  the  same  results  as  NUMIT.  It  is  written  in  Java.  See  [Davis,  et  al.,  2007]. 

DICTAT 

The  Spenvis  tool  was  developed  at  QinetiQ  (then  DERA).  See  on-line  documentation  and 
[Rodgers,  et  al.,  1998]. 

Anomaly  Investigations 
Violet  and  Frederickson 

An  anomaly  investigation  that  doesn't  get  the  attention  it  deserves  was  done  by  Mike  Violet  for 
Robb  Frederickson  [Violet  and  Frederickson,  1993]. 

SEAES-GEO 

This  expert  system  [O'Brien,  2009]  uses  the  12  hour  average  of  the  electron  flux  over  2  MeV 
measured  by  GOES  to  predict  the  anomaly  risk.  The  formula  is  based  on  the  correlation  of 
electrostatic  discharges  measured  by  SCATFIA  with  the  12  hour  average  of  the  over  1810  keV 
energy  channel  of  the  SC3  instrument.  The  tool  does  include  a  mapping  from  the  GOES  location 
to  the  spacecraft  location. 
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6.  MEMO  SUMMARIZING  RESOURCES  FOR  AURORAL  CHARGING 
COMPUTATIONS 


leidos 


Date:  February  4,  2014 

To:  Dale  Ferguson  and  Adrian  Wheelock 

From:  Victoria  Davis  and  Myron  Mandell 

Subject:  Summary  of  Resources  Regarding  Auroral  Charging  Environment  Characterization 
and  Charging  Computation 

From  a  surface  charging  perspective,  the  auroral  environment  is  poorly  characterized.  There 
exists  quality  data  that  could  be  used  to  fill  this  gap. 

DMSP 

The  best  correlation  study  between  environment  measures  and  spacecraft  potentials  we  are  aware 
of  is  that  done  by  Paul  Anderson1.  From  an  examination  of  12  years  of  DMSP  data,  he 
concluded  that  the  DMSP  spacecraft  charged  more  than  100  V  negative  under  the  following 
conditions. 

•  The  spacecraft  was  in  darkness; 

•  The  plasma  density  was  less  than  1 04  cm"3;  and 

•  There  was  a  high  integral  number  flux  (>  1 08  electrons  cm'2  s'1  sM )  of  high  energy  (>1 4 
keV)_electrons. 

Mengu  Cho  also  did  a  statistical  study  of  DMSP  environment  data.2  He  used  data  from  the 
environment  sensors  for  three  selected  years,  filtered  it  for  availability  of  full  spectra  and  a 
quality  measure,  and  fit  the  spectra  to  a  double  Maxwellian.  He  then  categorized  the  fits  by 
resulting  densities  and  currents  and  determined  the  frequency  of  the  electron  current  from  the  fit 
exceeding  the  ion  current  from  the  tit  as  a  function  of  year,  season,  flux  of  high  energy  electrons 
(based  on  fit),  and  pole.  With  additional  assumptions,  he  then  used  the  results  in  MUSCAT  to 
construct  a  frequency  of  charging  estimate.  He  appears  to  have  not  included  any  information 
regarding  the  measured  charging. 

David  Cooke  has  proposed  a  worst  case  environment  for  charging  calculations.  It  is  based  on  a 
fit  to  a  measured  environment  during  an  auroral  charging  event.  If  I  remember  correctly,  the  low 
energy  plasma  density  was  then  reduced  from  the  measured  value. 

Freja 

The  Freja  team  recorded  39  charging  events  at  high  altitude  over  the  northern  aurora.  The 
spacecraft  passed  over  the  southern  aurora  at  much  lower  altitude,  but  no  data  was  recorded. 1 

•  No  charging  events  were  observed  for  cold  plasma  density  above  2  x  1 03  cm"3. 

•  The  record  potential  was  nearly  -3  kV,  during  which  the  measured  ion  density  was  about 
30  cm"3. 

•  Most  events  occurred  when  the  spacecraft  was  sunlit.  However,  the  ionosphere  below 
was  in  darkness. 
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POLAR  calculations  using  a  representation  of  the  Freja  geometry  and  materials  and  fits  to 
measured  environments  were  done.4  The  calculations  predicted  much  less  charging  than 
observed.  One  known  issue  with  these  calculations  is  that  the  power  law  component  of  the 
Fontheim  was  extended  down  to  0.5  V,  well  below  the  energy  range  of  the  measured  spectrum. 
Thus  the  calculations  had  large  numbers  of  low  energy  electrons  that  were  (presumably)  not 
actually  present. 

Possible  Charging  Index  Studies 

The  Freja  calculations  could  be  redone  using  improved  fits  and  Nascap-2k. 

A  comparison  of  measured  chassis  potential  and  calculated  chassis  potential  of  DMSP  for  a 
selection  of  events  could  be  done,  with  the  goal  of  arriving  at  a  proposed  worst  case 
environment. 

A  study  similar  to  that  we  did  with  the  LANL  data  w  ith  the  goal  of  determining  the  best  spectral 
shape  to  use  for  charging  studies  would  be  very  interesting.  Such  a  study  will  be  most  efficient  if 
it  involved  Paul  Anderson.  The  last  time  such  a  study  was  contemplated,  it  was  determined  that  it 
would  take  significantly  more  than  the  funds  presently  allotted  to  us  for  the  charging  index  study. 
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7.  SEMI-ANALYTIC  AND  PIC  (PARTICLE-IN-CELL)  METHODS  FOR 
QUANTIFYING  CHARGING  IN  DENSE,  COLD  PLASMA 

The  following  paper  was  prepared  for  the  13th  Spacecraft  Charging  Technology  Conference  in 
Pasadena,  CA  in  June  2014. 

Semi-analytic  and  PIC  (Particle-in-cell) 

Methods  for  Quantifying  Charging 
in  Dense,  Cold  Plasma 

V.A.  Davis.  MJ,  Mandelk  D,L.  Cooke,  D.C  Ferguson 


Abstract— Researchers  use  several  techniques  for  the 
calculation  of  sheath  structure  unit  surface  currents  in  short 
Debye  length  plasma  for  the  purpose  of  calculating  surface 
potentials  of  complex  spacecraft.  These  techniques  include  the 
alignment  of  analytic  results  fur  simple  geometries,  semi- 
analytic  methods  that  use  Poisson's  solution  with  an  analytic 
charge  density  formula  in  conjunction  with  particle  tracking. 
hybrid-HC  (particle  in  cell)  techniques  in  which  the  charge 
density  is  determined  by  tracking  inns  u  ml  assuming  Boltzmann 
electron  densities,  and  full  PIC.  Hybrid- PIC  and  full  PfC 
Techniques  can  be  applied  to  buth  sialic  and  dynamic  plasmas. 

We  compare  the  various  technique*  available  in  the  plasma 
modeling  code  Nitsiurp-2k,  along  with  the  sinnlytic  approach  used 
in  the  analysis  tool  F  W B.  We  review  the  strengths,  weaknesses, 
and  limitations  of  the  available  models.  In  the  appropriate  limits, 
each  approach  gives  the  analytic  result — within  the  accuracy  of 
the  calculation.  W  hen  used  under  conditions  outside  the  limits  of 
the  approximations,  results  are  not  reliable  and  may  be 
misleading.  At  the  same  time,  the  PIC  and  hybrid  PIC 
approaches  can  be  misleading  when  phase  space  is  not 
adequately  sampled. 

When  used  within  their  range  of  applicability .  these 
techniques  arc  powerful  tools  in  assessing  charging  in  short 
Debye  length  plasma. 

A  e  y hj u nts — ch arfti ng;  tmitelhig;  tfcnw  ptasttui 

I.  CHAKCilNUMi  l  NANISMS 

This  paper  surveys  the  various  appmaches  used  to  model 
spacecraft  charging  in  cold,  dense  plasma,  such  as  found  in 
low- Earth-orbit.  The  range  of  plasma  properties  under 
consideration  are  listed  in  Tabic  Tin  these  plasmas,  with  a  few 
exceptions  (such  as  when  thin  dielectrics  lead  to  high  surface 
capacitances},  the  current  is  high  enough  that  spacecraft 
surface  potentials  adjust  to  changes  in  the  spacecraft  and 
environment  within  milliseconds.  Dynamic  effects  are 
associated  with  wake  effects,  result  from  switching  of  solar 
array  s  on  and  off  and  occur  in  dynamic  ex per i merits.  However, 
the  main  interest,  and  the  focus  of  this  paper,  is  the  calculation 
of  steady-state  potentials. 
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In  low  energy  plasma,  the  net  surface  current  is  dominated 
by  the  incident  charged  particles.  This  is  different  from  higher 
energy,  more  tenuous  plasmas,  such  as  geosynchronous 
substorm  and  auroral  environments.  In  cold  plasma  the  energy 
of  the  incident  charged  particles  is  low.  so  few  secondary 
electrons  arc  created,  and  (except  in  the  more  tenuous 
environments)  even  the  ion  t  henna  I  current  exceeds  the 
phoioemilted  electron  current. 

In  low -Earth  orbit  the  spacecraft  velocity  is  on  ihe  order  of 
7  km  sec.  slow  compared  with  the  electron  l henna l  speed,  but 
fast  compared  u  ftb  the  thermal  speed  of  an  atomic  oxygen  ion. 
The  ion  and  electron  densities  immediately  behind  the 
spacecraft  arc  dramatically  reduced  because  on  l>  tons  with 
velocities  comparable  to  the  spacecraft  velocity  can  exist  in 
lifts  deep  wake*  region.  The  electron  density  would  be  near 
uniform  were  it  not  for  the  negative  potential  due  to  electron 
space  charge  from  the  excess  electrons.  In  this  mesothermal 
plasma,  the  ion  current,  and  to  a  much  lesser  extent  the  electron 
current,  depends  on  the  angle  of  the  surface  w  ith  respect  to  the 
motion. 

Sign  i Hearn  charging  in  cold  plasma  is  driven  by  spacecraft- 
generated  voltages.  The  resulting  potentials  are  of  the  same 
order.  A  common  example  of  spacecraft -generated  voltages  is 
the  positively  biased  interconnect  and  bus  bars  of  the  solar 
arrays  that  are  exposed  on  most  spacecraft.  These  eiccirodes 
colled  electrons  from  Lhe  surrounding  plasma,  driving  the 
spacecraft  ground  negative,  generally  to  about  90%  of  the 
entire  solar-urray-generaied  voltage.  The  floating  potential  is 
deicnrftned  by  the  ion  current  to  other  exposed  conductive 
surfaces.  The  electron  cun  cut  collected  by  all  the  positive 
potential  conducting  surfaces  (such  as  the  interconnects)  is 
balanced  by  the  ion  current  collected  by  the  negative  potential 
conducting  surfaces. 

In  the  absence  of  spacecraft-generated  potentials  and  on 
insulating  surfaces,  lhe  surface  potential  is  determined  by  the 
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balance  between  the  attracted  species  and  the  repelled  species. 
Thus  surfaces  float  negative  a  few  times  the  plasma 
temperature,  so  that  the  attracted  ion  current  is  balanced  by  the 
repelled  electron  current.  Under  unusual  conditions,  potentials 
on  the  order  of  the  ion  kinetic  energy  can  develop  in  order  to 
repel  the  ion  flux, 

II.  A  VAI I  AH!  l  T ECHNIQl 

I  he  doselv  related  techniques  of  current  balance  and 
charging  are  used  to  compute  the  equilibrium  potential  of 
surfaces  in  cold,  dense  plasma  In  current  balance*  the  net 
current  as  a  function  of  surface  potentials  is  computed  in  order 
to  find  the  zero  of  the  current- volt  age  relation.  When  the 
charging  technique  is  used,  initial  surface  potentials  are 
assumed.  The  resulting  surface  currents  and  an  assumed 
capacitance  are  used  to  determine  the  change  in  surface 
potentials.  The  process  is  then  iterated  until  the  zero  current 
condition  is  readied.  In  both  approaches,  the  current  as  a 
function  of  surface  potential  is  the  central  parameter. 

The  current  computations  depend  on  how  the  plasma  is 
modeled.  There  are  lour  basic  approaches,  along  with  several 
variations. 

■  Analytic  currents 

0  Currents  computed  using  a  formula 

*  Particle  T rack  mg/ A  nalytic  space  charge 

c  Poisson's  equation  solved  with  space  charge  specified 
by  an  analytic  formula 

g  Current  of  one  or  both  species  computed  by  tracking 
macroparticles. 

*  Hybrid  PIC 

o  Poisson's  equation  solved  with  space  charge  given  by 
ion  densities  from  macropartide  tracking  and  electron 
densities  from  an  analytic  formula. 

;  Current  of  one  or  both  species  computed  by  tracking 
macroparticles 

*  PIC 

Each  of  these  approaches  has  benefits  and  drawbacks. 
While  analytic  approaches  arc  faster  and  more  stable,  the 
underlying  approximations  may  or  may  not  be  applicable,  Self 
consistent  PIC  calculations  include  all  the  physical  phenomena, 
but  require  care  to  adequately  represent  the  phase  space 
distribution  function  and  are  computationally  intensive. 

,T  Analytic  currents 

In  this  approach,  formulas  for  plasma  currents  are  derived 
by  assuming  an  idealized  geometry  and  plasma  and  calculating 
an  exact  answer.  A  major  advantage  is  that  once  a  formula  has 
been  developed,  applying  it  to  a  new  system  is  quick  and  easy. 
The  applicability  of  a  formula  to  other  conditions  is  limited  by 
how  close  the  system  under  consideration  is  to  the  ideal. 
Generally  these  formulas  depend  on  local  values  such  as  the 
surface  potential  and  orientation  with  respecl  to  ihe  direction  of 
motion,  and  sometimes  the  surface  electric  field.  Here  we 
discuss  two  common  and  one  not-so-common  approximations. 
In  the  plasmas  of  interest,  in  the  absence  of  spacecraft  driven 
changing  surface  potentials,  it  is  reasonable  to  assume  that  the 
plasma  is  at  equilibrium  and  therefore,  is  Maxwellian, 


In  the  absence  of  motion,  if  a  negative  potential  collecting 
surface  is  very  large  compared  to  the  Debye  length,  the 
attracted  current  is  the  ion  thermal  current  and  the  repelled 
current  Is  the  electron  thermal  current  reduced  by  the 
Boltzmann  factor  This  condition  can  be  used  to  determine  the 
fl  oat  i  ng  po ten  t  i  a  i  o  f  an  i  so  I  ated  obj  ect , 


for  an  oxygen  plasma,  the  solution  to  this  equation  is 
$  =  -5,  ]  4  T,  five  times  the  temperature. 

When  the  collecting  surface  is  in  motion,  ihe  thermal 
current  is  multiplied  by  a  velocity  dependent  angular  factor  on 
the  ram  side.  The  simplest  approximation  is  to  use  the  ion  ram 
current  times  a  cosine  of  the  angle  with  respect  to  11  ow  for  the 
ion  current  and  the  same  factor  as  above  for  the  electron 
current. 


for  a  surface  directly  facing  the  rum  at  a  typical  low-Earth- 
orbit  speed  of  7.500  m/s  in  a  0,1  cV  plasma,  the  surface  is  at 
-0,1 05  VT  twice  the  temperature, 

2)  Flowing  Maxwellian 

We  can  calculate  the  current  to  a  surface  in  motion  in  the 
orbit-limited  approximation  (flowing  Maxwellian).  The  orb h- 
limited  approximation  is  commonly  used  when  modeling 
tenuous  plasmas.  It  is  strictly  correct  for  a  l/r  potential  field 
with  no  physical  or  electrostatic  barriers. [  1 1  ft  is  appropriate  for 
cases  in  which  the  Debye  length  is  significantly  larger  than  the 
spacecraft  size,  and  therefore,  while  it  is  often  useful  as  a  first 
estimate,  the  results  must  be  confirmed  by  other  techniques 
when  used  in  dense,  cold  plasma.  The  current  is  given  by  the 
following  expression  [2j 


where  the  integration  variable  E  represents  the  energy  at  the 
surface,  the  integral  limit  L  is  0  for  the  repelled  species  and  |<jj| 

for  the  attracted  species,  the  fraction  in  parentheses  represents 
the  orbit  limited  enhancement  or  reduction  of  the  current,  the 
upper  sign  is  for  ions  and  the  lower  sign  is  for  electrons.  U  is 
the  velocity  of  the  plasma,  and  \j/  is  the  angle  between  the  flow 
vector  and  the  incident  velocity  at  infinity,  which  can  be  related 
to  the  velocity  at  which  the  particle  strikes  the  surface.  In  the 
a  bse  nee  o  f  m  ot  ion .  t  h  is  form  u  I  a  red  uces  t  o  a  Max  we  1 1  fa  n . 

In  the  orbit-limited  approximation,  the  potential  of  a  I  cm 
radius  conductive  stationary  sphere  in  a  lb11  m'\  0.1  eV 
oxygen  plasma  (Xq  =0.7  cm)  is  -0.36  V.  less  negative  than  in 
the  planar  approximation.  The  floating  potential  of  the  same 
sphere  moving  at  a  typical  low  Earth  orbit  velocity  of 
7.5  km /sec  is  -0,324  V.  slightly  less,  as  the  higher  ion  current 
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in  the  ram  more  than  makes  up  for  the  lower  ion  current  in  the 
wake. 


3)  Sheath 


In  dense  plasma  with  high  potentials,  the  disturbed  region 
is  referred  to  as  a  sheath.  In  the  sharp  sheath  edge 
approximation,  the  sheath  edge  represents  a  demarcation 
between  a  low  potential  exterior  region  containing  neutral 
undisturbed  plasma  and  a  high  potential  interior  region  from 
which  one  species  is  excluded.  The  charge  of  the  attracted 
species  in  the  sheath  region  balances  the  charge  on  the 
spacecraft  surfaces. 


Consider  a  simple  system  of  two  conductive  spheres  biased 
with  respect  to  each  other.  The  spheres  are  10  cm  in  radius  and 
their  centers  are  40  cm  apart.  Fig.  1  shows  the  potentials  in  a 
plane  through  their  centers  with  the  sphere  potentials  set  to 
+50  V  and  -50  V.  In  this  calculation,  the  plasma  is  10"  m"\ 
0.2  eV  giving  a  Debye  length  of  1.05  cm,  short  compared  with 
the  sphere  size  and  separation.  The  ions  are  atomic  oxygen  and 
the  spheres  are  stationary.  Within  each  sheath,  the  repelled 
species  is  excluded  leaving  behind  the  charge  of  the  attracted 
species  to  balance  the  charge  on  the  sphere.  The  overlapping 
sheaths  extend  about  10  cm  out  from  the  spheres.  Outside  of 
that  region  the  potentials  are  low. 

When  a  sheath  forms,  the  current  to  an  object  can  be 
estimated  as  the  current  to  its  sheath.  This  is  the  approach  used 
by  EWB  [3]  (Environment  WorkBench).  In  this  approximation, 
each  component  is  surrounded  by  a  sheath  of  radius  rsc  where 
rsc  is  the  solution  to  the  equation 


r  \ 

4/3 

(f  \3/4 

(  \ 

\ 

=  0.8356 

r  f\c 

_ 

a 

0.17 

j 

,U  J 

r 

\  SC  J 

/ 

4/3 


(4) 


where  a  is  the  radius  of  a  sphere  with  the  same  surface  area  as 
the  collecting  area  of  the  component.  This  formula  is  based  on 
earlier  work  by  Parker.[4]  The  sheath  surfaces  are  then  divided 


Fig.  1.  Illustration  of  potentials  on  plane  through  center  of  stationary 
two-sphere  object  with  spheres  are  +50  V  and  -50  V,  in  101 1  m  \  0.2 
eV  oxy  gen  plasma. 


into  several  parts.  Portions  of  the  sheath  that  overlap  are 
discarded.  The  current  to  each  sheath  segment  is  determined  by 
the  low  potential  formula  above.  This  approach  works  well 
when  the  problem  is  sheath  dominated.  In  this  model,  the 
floating  potentials  of  the  two  spheres  are  -19.6  V  and  +80.4  V. 

B.  Tracked  currents  in  potentials  computed  with  analytic 
space  charge 

A  more  computationally  intensive  approach  is  to  numerically 
solve  Poisson’s  equation  using  an  analytic  formula  for  the 
space  charge  and  then  track  charged  particles  in  the  resulting 
space  potentials.  One  such  approximation  that  is  applicable  at 
both  high  and  low  potentials  is  the  following  non-linear 
formula[5] 

_P_  = _ <|>  ma.r(l.C(<{>.  E)) 

e„  \D  l+V4n|<|)/7'|3/2 

C’(<t>,  E)  =  mm((R/rf  ,3.545|<t>/7j3/2  ) 

(R/rf  =2.29\EXjT\'2y^™ 

where  the  symbols  refer  to  the  local  potential,  <|>,  and  the  local 
electric  field,  E.  This  function  smoothly  interpolates  between 
linear  Debye  screening  at  low  potentials  and  the  charge  density 
of  a  single  accelerated  and  converging  species  at  high 
potentials.  The  quantity  C,  which  is  a  function  of  the  local 
potential  and  electric  field,  accounts  for  the  increase  in  charge 
density  as  charged  particles  from  a  large  area  are  attracted  to  a 
small  region.  The  convergence  formula  was  developed  to  fit 
the  results  of  Langmuir  and  Blodgett  [6]  for  current  collection 
by  a  sphere.  In  a  dense  plasma,  when  the  spacecraft  velocity 
and  Earth’s  magnetic  field  have  minimal  effect  on  the  charge 
density  within  the  sheath,  this  formula  is  appropriate. 

Once  Poisson’s  equation  has  been  solved  for  space 
potentials,  surface  currents  can  be  determined  by  tracking 
macroparticles  either  from  a  sheath  edge  or  from  the  boundary 
of  the  computational  space  as  appropriate. 

When  a  sheath  exists,  it  can  be  more  convenient  to  track 
macroparticles  from  the  sheath  edge.  Macroparticles  can  be 
created  at  the  sheath  edge  and  then  tracked  in  the  computed 
potentials  and  any  specified  magnetic  fields  to  determine  the 
current  to  each  surface.  In  the  computation,  the  sheath  edge  is  a 
surface  at  a  specified  potential,  for  example,  <|>  =  ±T  In  2  .  The 
justification  for  this  choice  is  that  because  the  attracted  species 
is  absorbed  by  the  sheath,  only  the  inward  moving  component 
is  present,  comprising  half  the  ambient  density.  The  repelled 
species,  whose  density  satisfies  n(ty)  =  n  exp(-\<b/T\),  must  also 

be  at  half  the  ambient  density,  leading  immediately  to  this 
sheath  potential.  The  current  for  each  macroparticle  is  the 
thermal  current  for  the  sheath  area  it  represents.  The  simplest 
assumption  for  the  initial  velocity  of  the  macroparticles  is  the 
average  velocity  of  charged  particles  crossing  the  sheath  edge. 


The  approximations  made  when  tracking  from  a  sheath  are 
not  valid  when  either  the  thermal  distribution  of  particle 
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velocities  or  the  spacecraft  velocity  is  important.  In  these 
situations  macroparticles  are  tracked  from  the  boundary  of  the 
computational  grid.  Enough  macroparticles  need  to  be  tracked 
to  capture  the  distribution  function.  The  velocity  and  current 
carried  by  each  macroparticle  can  be  assigned  stochastically  or 
a  specified  portion  of  the  distribution  function  can  be  explicitly 
specified.  Macroparticles  can  be  distributed  so  that  they  have 
approximately  the  same  weight  or  macroparticles  with  less 
current  can  be  used  to  represent  the  edges  of  the  distribution 
function. 

A  geometric  wake  factor  can  be  added  to  the  charge  density 
and  sheath  current  formulas  to  account  for  spacecraft  (and/or 
plasma)  motion. 

The  stationary  two-sphere  example  discussed  above  can  be 
used  to  illustrate  this  type  of  model.  The  potentials  shown  in 
Fig.  2  were  computed  by  Nascap-2k  using  the  charge  density 
formula  of  (5).  The  current  to  each  sphere  was  computed  both 
by  tracking  macroparticles  from  the  sheath  edge  and  by 
tracking  particles  from  the  boundary  of  the  computational 
space  and  the  sphere  potentials  were  adjusted  to  find  the 
potentials  for  which  the  net  current  is  zero.  The  surface 
potentials  of  +0.25  V  and  -99.75  V,  shown  in  Fig.  2,  are  those 
that  result  when  tracking  from  the  sheath.  This  result  is  20  V 
different  than  that  obtained  for  the  same  sphere  size  and  plasma 
properties  using  the  analytic  model  of  sheath  current  above. 

When  tracking  from  the  sheath  edge,  the  current  to  each 
sphere  includes  only  the  attracted  species.  When  tracking 
macroparticles  of  both  species  from  the  boundary  of  the 
computational  space,  it  is  possible  to  find  a  floating  potential 
with  both  spheres  negative.  For  this  geometry,  the  floating 
potentials  of  -100.2  V  and  -0.2  V  were  found.  The  -0.2  V 
sphere  collects  both  ions  and  the  high  energy  tail  of  the 
electron  phase  space  distribution  function,  which  when 
combined  with  the  ion  collection  of  the  -100.2  V  sphere  gives 
zero  net  current. 

As  shown  in  Fig.  3,  at  lower  density,  the  large  (negative) 


TABLE  II.  TWO-SPHERF  FLOATING  POTENTIALS  (V)  COMPUTED 
USING  DIFFERENT  TECHNIQUES 


Debve  length  (cm) 

1.05 

3.32 

10.5 

Analytic  currents.  Sheath 

-90.9 

+9.1 

-19.6 

+80.4 

Analytic  currents,  Maxwellian 

-99.6 

+0.4 

-99.6 

+0.4 

-99.6 

+0.4 

Tracked  currents  from  sheath 

-99.75 

+0.25 

-98.1 

+  1.9 

-86.1 

+  13.9 

Tracked  currents  from  boundary 

-100.2 

-0.2 

-98.7 

+  1.3 

-88 

+  12 

sheath  gets  larger,  enveloping  the  positive  sphere.  In  order  to 
achieve  current  balance,  the  positive  sphere  is  at  a  higher 
potential,  here  13.9  V,  so  that  the  positive  potential  peeks 
through  and  electrons  can  be  collected.  For  this  simple 
geometry,  the  ratio  of  the  exposed  sheath  areas  matches  the 
ratio  of  plasma  thermal  currents. 

The  floating  potentials  for  the  two-sphere  example  for 
1011  m'\  1010  nr  and  109  m'\  0.2  eV,  oxygen  plasmas  are 
shown  in  Table  II.  The  results  using  the  Analytic  Currents, 
Sheath  model  matches  the  other  results  best  in  the  shortest 
Debye  length  case,  for  which  it  is  expected  to  be  most  accurate. 
For  this  example,  the  results  are  not  sensitive  to  whether  the 
currents  are  tracked  from  the  sheath  edge  or  the  boundary. 
When  the  spacecraft  is  moving  and  when  the  currents 
themselves  are  important  (rather  than  just  their  relative  values), 
it  can  be  important  to  track  the  macroparticles  from  the 
boundary. 

C.  Hybrid  PIC 

When  the  spacecraft  geometry  is  complex  enough  that 
analytic  charge  densities  are  inadequate,  the  charge  density  can 
be  given  by  the  sum  of  an  ion  density  determined  by  tracking 
macroparticles  and  an  electron  density  given  by  a  Boltzmann 
function.  The  equilibrium  electron  density  is  Boltzmann  for 
negative  potentials  in  the  absence  of  a  potential  barrier  and 
nearly  Boltzmann  for  small  barriers. 

As  the  ion  macroparticles  are  used  to  give  the  ion  charge 
density,  they  must  be  tracked  from  the  computational  boundary. 


Fig.  2.  Potentials  on  plane  through  center  of  stationary  two-sphere 
object  with  spheres  at  their  floating  potentials,  of  +0.25  V  and  - 
99.75  V,  computed  by  tracking  macroparticles  from  the  boundary  of 
the  computational  space,  in  10"  m*\  0.2  eV  plasma. 
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To  get  self-consistent  potentials  and  ion  densities,  there  are  two 
approaches.  Hie  ions  can  be  tracked  throughout  the  grid  and 
then  the  resultant  potentials  computed  and  the  process  iterated, 
possibly  including  sharing  of  ion  densities  between  iterations. 
Alternatively,  a  more  traditional  PIC  approach  can  lie  used  in 
which  the  computational  space  is  tilled  with  ion 
macro  part  ides.  The  maeroparticles  are  tracked  fora  short  time, 
the  potentials  are  computed,  additional  ions  are  injected  from 
the  computational  boundary,  and  the  process  repeated  until  ail 
transients  disappear.  In  both  eases,  the  time  required  to  reach 
equilibrium  depends  on  how  close  the  initial  potentials  are  to 
the  potentials  for  zero  Current.  As  mentioned  above  for  electron 
tracking  to  obtain  surface  currents,  the  initial  velocities  of  the 
macroparticles  used  to  represent  the  velocity  distil  but  ton  of 
phase  space  can  be  assigned  stochastically  or  a  specified 
petition  of  the  distribution  function  can  be  explicitly  specified. 
The  accuracy  of  the  potential  solution  land  therefore  the 
resulting  currents  and  zero  current  surface  potentials)  is  limited 
by  the  adequacy  of  sampling  of  the  ion  phase  space  distribution 
function. 

The  tracked  ion  macroparticles  can  be  used  to  compute  the 
surface  currents  as  well  as  ihe  spatial  charge  density.  The 
electron  surface  current  can  be  computed  using  any  of  the 
analy  tic  or  tracked  approaches  described  above. 

The  importance  of  velocity  resolution  can  be  illustrated  by 
considering  the  floating  potential  of  an  isolated,  conducting 
sphere.  While  an  isolated  sphere  or  cube  is  a  favorite  for  testing 
charging  codes,  its  symmetry  and  low  potential  make  it 
numerically  challenging.  We  nse  a  stationary,  isolated  I  cm 
radius  conductive  sphere.  In  a  tOlft  m  '  to  Iff'  m1.  0.1  eV, 
oxygen  plasma,  the  Debye  length  of  0,235  cm  to  2.35  cm  is 
comparable  to  the  sphere  size.  Above  we  obtained  a  floating 
potential  of -0.5!  V  in  (he  planar  limit  and  -0.36  V  in  the 
orbit- limited  current  collection  limit.  Table  III  gives  the 
floating  potential  of  this  sphere  computed  using  an  analytic 
charge  density  and  Hybrid  PIC  techniques  with  various 
representations  of  the  velocity  distribution. 

First,  we  use  the  analytic  charge  density  function  when 
computing  Ihe  space  potentials,  track  a  distribution  of  ion 
macroparticles  in  those  potentials  from  the  boundary  of  the 
computational  space,  and  use  the  Boltzmann  relation  to  specify 
the  electron  current.  When  the  Debye  length  is  larger  than  the 
sphere,  as  expected,  the  orbit  limited  result  is  obtained,  In  a 
shorter  Debye  length  plasma,  the  surface  potential  is  shielded 
by  the  attracted  ions,  so  that  ion  collection  is  less  than  orbit - 
limited,  and  the  sphere  floats  more  negative.  Fig.  4  shows  the 
potential  in  a  plane  through  the  center  of  the  sphere. 

If  the  same  calculation  ri  performed  using  the  Hybrid  PIC 
approach  with  a  minimal  distribution  of  initial  velocity  values 
of  the  ion  maeroparticles,  the  resulting  sphere  potential  is 
slightly  more  negative  and  the  space  potentials  seen  in  Fig.  5 
result.  In  this  calculation,  the  thermal  velocity  distribution  is 
represented  by  eight  evenly  weighted  maeroparticles  with  the 
average  thermal  velocity  plus  a  component  to  represent  ihe 
thermal  distribution  about  a  randomly  generated  coordinate 
axis.  That  the  lowest  potential  contour  is  somewhat  square,  like 


TABU-:  It! .  Fl BATING  POTENTIAI  ( V}  OP  SPHFIlli  COMH  ft  I  D  1  StNt  i 

UJPhtRtN  I  r LL'J  IN  WJL  Ir.S 


Debye  length  (cm) 

0  2 A 

0,74 

235 

Analytic  charge  density,  tracked  ion  current  and 
analy  lie  electron  current 

-0.4(J 

-0,37 

-0.36 

Hybrid  PTC  wrlh  poor  ion  velocity  resolution 

-0.40 

Hybrid  PIC  w  ills  heller  ion  velocity  restilulion 

-G.3B 

-0.37 

-0  36 

Hybrid  PIC  and  irnek  electrons  vvilh  poor  \  duchy 
resolution  for  current* 

-0.33  ‘ 

hQJO 

-0.28 

1  lybrid  EMC  and  track  electrons  with  better 
velocity  resolution 

■0  37 

-034 

-0.33 

the  computation  grid,  rather  than  circular,  like  the  sphere, 
suggests  that  the  results  may  not  be  correct.  Adding  additional 
ion  macroparticles  with  a  wider  range  of  initial  velocities, 
shown  in  Fig  6  gives  a  result  similar  to  the  analytic  charge 
density  calculation  (which  for  this  calculation  is  expected  to 
give  a  highly  reliable  result).  A  total  of  64  maeroparticles  were 
created  at  each  Initial  location  In  each  of  the  three  spatial 
directions,  these  maeroparticles  represent  the  10%  of  ions 
moving  fastest  in  the  negative  direction,  the  other  40%  moving 
in  the  negative  direction.  10%  of  the  ions  moving  fastest  in  the 
positive  direction,  and  the  other  40%  moving  in  the  positive 
direction. 

When  the  electron  current  is  computed  using 
maeroparticles,  the  adequacy  of  the  representation  of  the  high 
energy  portion  of  phase  space  becomes  very  important,  as  can 
be  seen  by  comparing  the  floating  potentials  computed  with 
poor  and  better  resolution.  In  the  first  of  the  calculations  in 
which  the  electron  current  is  determined  by  tracking  electron 
maeroparticles.  the  electron  thermal  distribution  is  represented 
in  the  same  manner  as  the  ion  distribution.  This  is  inadequate 
as  only  a  small  fraction  of  the  electrons  have  enough  energy  to 
reach  the  sphere.  In  the  better  resolution  calculation,  a  similar 
scheme  for  representing  the  distribution  is  used.  At  each 
emitting  point  on  the  boundary,  1,000  maeroparticles  are 
created.  The  electrons  distribution  is  divided  into  those  with  the 
2.5%  most  negative  velocities,  those  with  the  next  2.5%  most 
negative,  those  with  the  next  5%  most  negative,  those  w  ith  the 
next  5%  most  negative,  the  remaining  35%  with  negative 
velocities,  and  the  same  for  positive  velocities. 

a  PIC 

The  first  principles  approach  is  to  use  macro  particle 
tracking  of  all  species  to  determine  the  charge  density  for 
Poisson's  equation  as  well  as  to  determine  the  surface  currents. 
This  requires  even  more  computer  time  and  a  great  deal  of  skill 
to  avoid  numeric  pathologies.  There  is  an  extensive  literature 
of  techniques  to  help  speed  solutions  and  improve  stability 
developed  by  researchers  modelling  dynamic  plasma 
conditions.  The  calculation  is  inherently  dynamic,  so  averaging 
and  smoothing  are  generally  needed  to  determine  the  steady- 
state  solution. 

This  approach  becomes  necessary  when  investigating  novel 
effects  and  exploring  dynamic  conditions. 
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Potentials 


Fig.  4.  Potentials  in  a  plane  through  the  center  of  a  sphere  computed  using 
the  analytic  space  charge  density  of  (5)  at  the  floating  potential  of 
-0.37  V  in  a  101 1  m  \  0. 1  eV  plasma. 


Fig.  5.  Potentials  in  a  plane  through  the  center  of  a  sphere  computed  using 
the  Hybrid  PIC  method  with  minimal  velocity  resolution  at  the 
floating  potential  of -0.40  V  in  a  1011  m'\  0. 1  eV  plasma. 


Fig.  6.  Potentials  in  a  plane  through  the  center  of  a  sphere  computed  using 
the  Hybrid  PIC  method  with  better  velocity  resolution  at  the  floating 
potential  of  -0.40  V  in  a  1011  m'\  0.1  eV  plasma. 


III.  Conclusions  and  Implications 

We  have  surveyed  a  selection  of  methods  to  compute 
currents  when  solving  for  floating  potentials  of  surfaces  in 
dense,  cold  plasma  such  as  observed  in  low  Earth  orbit.  The 
methods  range  from  analytic  to  full  PIC.  Analytical  methods 
are  easy  to  apply  and  give  results  quickly.  However,  they  are 
only  applicable  when  the  underlying  approximations  are  valid. 
At  the  other  extreme,  full  PIC  calculations  are  applicable  to 
complex  geometry  and  dynamic  conditions.  However,  they 
require  large  investments  of  computational  resources  and  a 
great  deal  of  skill  to  obtain  reliable  results. 

In  practice,  computations  generally  combine  analytic 
methods  with  PIC  techniques.  With  an  understanding  of  the 
underlying  approximations,  reliable  results  can  be  obtained 
with  the  minimum  necessary  resources. 

The  extent  to  which  the  approximations  are  valid  is  highly 
dependent  on  the  system  and  the  question  at  hand.  To  gain 
confidence  in  the  results,  it  is  best  to  use  the  simplest 
applicable  models  to  gain  insight,  and  then  use  more  complex 
approaches  to  verily  the  applicability  of  the  approximation. 
When  calculations  are  dependent  on  adequate  resolution,  they 
can  be  checked  by  verifying  that  the  results  are  the  same  with 
improved  resolution. 
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